/*
 *  Matrix.cpp
 *  T3nsors
 *
 *  Created by Michael Barriault on 09-10-27.
 *  Copyright 2010 University of Guelph. All rights reserved.
 *
 */

#include "Matrix.h"
#include "Macros.h"
#include "Others.h"
#include <cstdarg>

Matrix::Matrix(void) {
	return;
}

Matrix::Matrix(Domain iO, int iN, std::string iid) {
	O = iO;
	N = iN;
	id = iid;
	T.assign(N, Vector(O,N,id));
	FOR(n,N)
		T[n].id = T[n].id + itos(n);
}

Matrix::~Matrix(void) {
	T.clear();
}

Vector& Matrix::operator()(int i) {
	return T[i];
}

Scalar& Matrix::operator()(int i, int j) {
	return T[i](j);
}

double& Matrix::operator()(int m, int n, int a, ...) {
	va_list va;
	va_start(va, a);
	int t=0;
	FOR(i,O.N) {
		int A = a;
		for ( int j=i-1; j>=0; j-- )
			A *= O[j];
		t += A;
		a = va_arg(va, int);
	}
	va_end(va);
	return T[m](n)[t];
}

Matrix diff(Vector x) {
	Matrix dx(x.O, x.O.N, (std::string)"d"+x.id);
	FOR(p,dx.N)
		dx(p) = diff(x(p));
	return dx;
}

Scalar div(Vector x) {
	Matrix Dx = diff(x);
	Scalar dx(x.O, (std::string)"div"+x.id);
	FOR(o,x.O()) FOR(a,x.N) dx[o] += Dx(a,a)[o];
	return dx;
}

Vector curl(Vector x) {
	Matrix Dx = diff(x);
	Vector dx = Vector(x.O, x.N, (std::string)"curl"+x.id);
	dx(0) = Dx(2,1)-Dx(1,2);
	dx(1) = Dx(0,2)-Dx(2,0);
	dx(2) = Dx(1,0)-Dx(0,1);
	return dx;
}

